SCIENTIFIC 

REPORTS 




The Cell Cycle Switch Computes 



SUBJECT AREAS: 

COMPUTATIONAL 
BIOLOGY 

INFORMATION THEORY AND 
COMPUTATION 

SYSTEMS BIOLOGY 

CELL SIGNALLING 



Received 
6 July 2012 

Accepted 
29 August 2012 

Published 
1 3 September 201 2 



Approximate Majority 

Luca Cardelli' & Attila Csikasz-Nagy 2 



Microsoft Research, 7 J J Thomson Ave Cambridge CB3 OFB, UK, 2 The Mic rosoft Research-University of Trento Centre for 
Computational and Systems Biology, Piazza Manufattura, Rovereto-38068, Italy, 3 Randall Division of Cell and Molecular 
Biophysics, and Institute for Mathematical and Molecular Biomedicine, King's College London, London, SE1 1 UL, UK. 

Both computational and biological systems have to make decisions about switching from one state to 
another. The 'Approximate Majority' computational algorithm provides the asymptotically fastest way to 
reach a common decision by all members of a population between two possible outcomes, where the decision 
approximately matches the initial relative majority. The network that regulates the mitotic entry of the 
cell-cycle in eukaryotes also makes a decision before it induces early mitotic processes. Here we show that the 
switch from inactive to active forms of the mitosis promoting Cyclin Dependent Kinases is driven by a 
system that is related to both the structure and the dynamics of the Approximate Majority computation. We 
investigate the behavior of these two switches by deterministic, stochastic and probabilistic methods and 
show that the steady states and temporal dynamics of the two systems are similar and they are exchangeable 
as components of oscillatory networks. 
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At the core of the biochemical networks that regulates the cell cycle there is a switch that triggers an 
irreversible transition from one critical stage to the next 1 . To understand the general properties of this 
switch across different species one may simplify the molecular interactions found in real organisms and 
abstract over the molecular species that are used by different organisms for the same functions 2,3 . In that abstract 
sense, the cell cycle transition to M phase has been shown to employ a universal switching mechanisms in all 
eukaryotes 4,5 . The mitotic entry regulator Cyclin Dependent Kinase (Cdk) activates its own activator Cdc25 and 
inhibits its inhibitor Weel and with these two positive feedback loops the system can abruptly switch from low to 
high Cdk activity state. The structure and properties of this molecular switch have been widely studied 6 8 , and 
successful mathematical models of the switch and its functional context (the cell cycle oscillator) have been 
built 9,10 . There is evidence that related networks are used in other cell cycle transitions 11,12 . 

From the point of view of dynamical systems, the basic features of the cell cycle switch are essentially 
understood: non-linear kinetics and positive feedbacks induce bistability (needed for switching) and result in 
hysteresis (needed to resist switch-back) 13 ~ 15 . Such kinetic behavior, however, could be achieved by many different 
networks, while in nature we find a specific universal structure of the network. For example, the required positive 
feedbacks are produced by double-negative and double-positive loops together 16,17 . What is special about this 
structure of the universal network, apart from its required function? How does the system handle noise and how 
does it relate to other non-biological switches? 

To answer these questions we shift our perspective from dynamical systems to computing systems. We ask: 
what does the cell cycle switch compute, and how does it compute it? We find that a fundamental algorithm from 
distributed computing 18,19 , the Approximate Majority (AM) algorithm 20 matches both the networks structure and 
the kinetics of the cell cycle (CC) switch; therefore, as we shall describe, the network structure derives from the 
algorithm, subject to biological constraints. The AM algorithm computes the majority of two finite populations by 
converting the minority population into the majority population, so that a single population results; it uses a third 
'undecided' state of the population, from where autocatalysis can drive the individuals into either of the final 
states. The algorithm has been studied in the context of Population Protocols 21 : a model of computation and 
related algorithms for resource-limited interacting agents (such as networks of environmental sensors). It not 
only switches a majority into a totality, but it does so in a particular way that is (1) fast: with high probability only 
0(N log N) binary interactions (consult the Methods section for asymptotic notation) are required for a popu- 
lation of size N ( 20 Theorem 1), (2) reliable: above a certain threshold the true initial majority wins with high 
probability ( 20 Theorem 2), and (3) robust: the algorithms still functions correctly in presence of a subpopulation 
of the order of sqrt(N) that behaves incorrectly or even 'deviously' ( 20 Theorem 4). Moreover, the algorithm is 
asymptotically optimal in the number of reactions required to switch a majority into a totality, because at least 
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order of N log N interactions are required for each element of the 
population to interact, directly or indirectly, with every other mem- 
ber. An interesting question is how close nature can come to this 
optimum. 

The AM algorithm can be described chemically in terms of just 4 
catalytic and autocatalytic reactions (equation 1). The comparison 
between AM and CC, however, is not straightforward. To carry it out 
we need to extend our understanding of both systems: in the com- 
puting literature it is not common to consider AM in terms of 
dynamical systems, which are usually continuous, and in the bio- 
logical literature is not common to consider CC in terms of com- 
putational systems, which are usually discrete. The AM network 
structure is different from (simpler than) the CC structure, but we 
argue that the differences are dictated by biological constraints, and 
that both the computational and the dynamical properties of AM are 
preserved by the transformations required to turn it into CC so that 
the constraints are satisfied. In particular, we show that (1) the struc- 
ture of AM well implements an input-driven switching function (in 



addition to the known majority function), (2) the structure of CC 
well implements a input-less majority function (in addition to the 
known switching function), (3) the structures of AM and CC are 
related, and an intermediate network shares the properties of both, 
(4) the behaviors of AM and CC in isolation are related, (5) the 
behaviors of AM and CC in oscillator contexts are related, (6) a 
refinement of the core CC network, known to occur in nature, 
improves switching performance and brings it in line with AM per- 
formance. 

Results 

We investigate and compare four networks that perform switching 
between a population x and a population}'. In all our systems (Fig. 1) 
x and y are converted to each other by reactions that are directly or 
indirectly driven by x andy, in such a way that each system contains 
two positive feedback loops that may also act through intermediaries 
or by double inhibition. 
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Figure 1 | Switch kinetics in absence of external control. We study the convergence time (time to reach a stable state) of switching networks; in each 
column, one of two stable states can be reached from the initial conditions. All the reaction rates are set equal and all the species start at equal quantities 
(details are provided in Supplementary Materials). Columns A-D concern the DC, AM, SC, CC switches respectively. Row 1 gives a depiction of the 
networks from which one can precisely recover the chemical reaction networks according to our notation. A catalytic reaction is represented by a circle on 
top of an arrow; for example, the top right of (F) contains the reaction b+z-^z+y with catalyst z. A reaction like x+y—>y+y is said to be autocatalytic. The 
full diagram (F) depicts two catalysts, z and r, acting on species x and y through a shared intermediary b, representing the 4 reactions x+ z-^>z+ b, 
b+z-^z+y, y+r^r+b, and &+r— »r+x We use a more compact pinched-arrow graphical notation (E) to represent the same network as (F), hiding the 
intermediary species b that is assumed not to enter any other reaction. Note that if multiple catalysts act on the same pinched arrow (as in Fig. 2), they all 
act on the hidden intermediary in the same way. Row 2 contains deterministic simulations for the mass action ODEs of the respective systems for four 
values of the initial discrepancy between xand y, from 10% to 0.01% (not meaningful for DC as the system would rest at its initial setting), (Y-axis scale on 
the right) . Row 3 shows sample stochastic simulations consisting of individual traces (black lines) of the Gillespie algorithm. The background heat maps, 
shown in logarithmic scale, give the probability Pr(x, -\ tj) that a system will have x, molecules ofxat time tj. The horizontal axis is time, and the vertical axis 
is concentration (row 2) or number of molecules (row 3). Subscripts V are for stochastic variables and 'p' for probabilistic variables. 
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Figure 2 | Switch equilibrium in presence of external control. We study the steady state response of our switching networks. Row 1 shows the four systems 
of figure 1 extended with external controls sx (switch-to-x) and sy (switch-to-y) in gray. On row 2, the horizontal axis represents the value of the sx input, 
and the vertical axis the value of the resulting x output at steady state. The value of sy remains fixed throughout; as in figure 1 all other initial values and 
rates are set equal. In each plot, the black line is the bifurcation diagram for the mass action ODEs of the systems (solid = stable, dashed = unstable steady 
states), showing that each system except the first one exhibits hysteresis; the axes represent concentrations with values that match the ovelayed plots 
(arbitrary units). The heat map in the background details the discrete probability distribution (in log scale) of a system composed of 5 molecules of each 
species, except for sy= 3 and with sx varying in 0.. 15. A point (sx h Xj, z*) in the heat map, for 0 S 1, means that at equilibrium (after a sufficiently long 
time) the discrete probability of the system being found with Xj molecules of species x on an input of sx^ molecules of species sx is equal to z*. The noisy red 
line is a single run of a steady-state stochastic simulation with maximum x= 1 50 and with sy= 30; sx is increased slowly to obtain the lower trajectory and 
the jump up, and is then decreased slowly to obtain the upper trajectory and the jump down. 



We consider first a system derived from AM by removing the 
intermediate species b, leading to two coupled positive feedbacks that 
works by Direct Competition (DC) between x andy. This network is 
depicted in figure 1A.1 and corresponds to the two autocatalytic 
reactions x + y — > y + y and y + x — » x + x. (A black ball over a 
reaction arrow represents presence of a catalyst.) The AM 
(Approximate Majority) network is depicted by an abbreviated gra- 
phical notation in figure 1B.1, which is explained in figure 1E,F. 
Namely, a pair of 'pinched arrows' between a pair of species such 
as x,y represents four reactions with a common intermediate spe- 
cies b (specific to the pair x,y): the intermediate species is omitted. 
The AM network therefore consists of 4 reactions (of mass action 
kinetics): 

x+y -*y + b catalytic in y 
b+y^y+y autocatalytic in y 
y + x^>x+b catalytic in x 
b + x— >x + x autocatalytic in x 

for some intermediate b. The pinched-arrow motif is suggestive of 
multi-site modification (e.g.: phosphorylation/dephosphorylation) 22,23 . 

In nature it is not common to find direct autocatalytic reactions, 
and especially mutually autocatalytic reactions like in AM. We can 
however replace them with mutually catalytic reactions that are not 
directly autocatalytic. We introduce intermediate catalysts r, acti- 
vated by x (from a reservoir p) and activating x, and z, activated by 
y (from a reservoir w) and activating y. This transformation results in 
the network SC (Simply Catalytic) of figure 1C.1, where we use 
pinched-arrow transitions for the new species too (other choices 
might be suitable). Here we need to introduce additional catalysts 
f,s to give a threshold for the self-activation reactions of x and y 
respectively. 



Natural systems that switch a molecule between two states will 
normally have that molecule being active in one state and inactive in 
the other. However in AM, and still in SC, both states x and y are 
active and catalyze different reactions. Bifunctional enzymes exists in 
various organisms 24,25 , but in the case of the cell cycle switch there is 
no evidence that Cdk would have separate type of activities in its 
different forms. To move closer to CC we further transform SC in 
figure 1D.1: we remove the catalytic activity from y, and we are 
therefore left with z activating y, which is inactive. To deactivate y 
we must therefore deactivate z; this function is taken over by x, which 
on one hand activates itself through r, and on the other hand deac- 
tivates y through z (note that the role of s also has reverted compared 
to its role in SC). This final transformation of the network is certainly 
non-obvious, but the network so obtained turns out to be the char- 
acteristic network of the cell cycle (CC) switch 4,8 , with x and y cor- 
responding to active and inactive forms of cyclin dependent kinases, 
and z and r resembling active forms of Weel and Cdc25 respectively 
(see Supplementary Materials for the exact correspondence). 

AM and CC as isolated systems. We now compare the performance 
of these four switches by deterministic 26 , stochastic 27 , and proba- 
bilistic analysis 28 . We first examine the convergence time of the 
switches: the time to switch from an unstable initial state to one of 
two stable states. In our analysis, all the rates are set to 1.0 and the 
populations start at similar levels. The second row of Figure 1 
contains simulations of the deterministic ODEs of the systems, for 
different initial small, but non-zero, differences between x andy. The 
black traces in the third row of Figure 1 are individual stochastic 
simulations for systems with molecular counts of the order of 
thousands of molecules, starting with all species at the same level 
including x=y. The deterministic system would not show the break 
of the symmetry from a totally balanced initial condition, while 
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stochastic noise is enough to move the systems out of this unstable 
steady state. The heat maps in the third row of Figure 1 are 
probability distributions for systems with low molecular counts, of 
the order of tens of molecules, again starting with all species at the 
same level. The three analysis techniques address different aspects 
and constraints. Deterministic simulation gives a characterization of 
the mean behavior of a system, but does not address noise charac- 
teristics. Noise can be readily observed in the stochastic simulations, 
but at very low counts and high noise individual traces give little 
information. The probability distributions, on the other hand, give 
a precise picture of a system at low count, but become computa- 
tionally unfeasible at higher count. 

Timing of transitions. The analysis of AM, SC, and CC indicates 
rapid switching and convergence to a stable state after an initial 
hesitation due to breaking of the symmetry of the initial condi- 
tions. If the initial conditions are not so symmetrical, the switch- 
ing happens more quickly: we essentially investigate the worst-case 
scenarios. We now discuss the performance of the individual 
switches. 

The Direct Competition switch (DC) (Fig. 1A.3) can rest in two 
different states, if all the molecules are ever found in state x or state y 
then the system remains in that state, because at least a single mole- 
cule from the other species is needed to catalyze reactions that move 
the system away from these states. However, the path to reach such a 
steady state from arbitrary starting conditions is a random walk, 
taking expected O(AP) molecular interactions to achieve, with each 
interaction moving the system randomly one way or the other. 
Moreover, any small perturbation of such a steady state has a chance 
to move the system to the other state, again by a random walk: the 
steady states are not robust to noise, thus the system is not bistable 
from a dynamical point of view. The deterministic ODE model of this 
system has all the right hand sides of the equations equal to zero, thus 
the deterministic system remains at its initial condition (not shown). 

The Approximate Majority switch (AM) is bistable (Fig. 1B.2-3), 
and the expected time to reach a stable state from any starting con- 
dition is with high probability less than 0(N log AO molecular inter- 
actions and 0(log AO time steps, and in fact approximately 3 log N in 
the worst case of equal starting conditions x = y 20 . The stable states 
are completely stable: once reached, no further reactions can happen. 
Moreover, if the initial state has a majority of either x or y exceeding 
w(sqrt(N) • logN) (a low percentage at high molecule numbers), then 
with high probability the system will settle in the state that has the 
initial majority 20 . This implies that the vicinity of the two stable 
states, where either x or y have a wide or total majority, is resistant 
to perturbations. 

The Simply Catalytic switch (SC) exhibits the same convergence 
behavior as AM modulo a small time factor (Fig. 1C.2-3). Note 
however that SC no longer stabilizes at the maximal level completely. 
The stable states are overwhelmingly but not completely stable 
because of the continued bias of the reverse reaction catalysts (s 
and i) which work against x. 

The Cell Cycle switch (CC) preserves much of the AM-like beha- 
vior, such as fast decision (Fig. 1D.2-3). In addition to the lack of 
complete stability, as in (SC), the switch is no longer fully symmet- 
rical: when y is in majority x approaches zero and y (not shown) 
reaches maximum, but when x is in majority it is counteracted by the 
fixed bias s that keeps some z active and constantly drives some x to y, 
thus x settles well below the maximum level. SC and CC have a bit 
slower transition time as well (approximately double the time, com- 
pared to AM), which might be caused by the same effect. These 
discrepancies between AM and CC are resolved in biological systems: 
in Discussion we examine extra feedbacks that appear in cell cycle 
systems 29 and maximize x levels. Possibly such feedbacks evolved 
exactly to allow a full activation of x, making CC dynamics even 
closer to that of the efficient system of AM. 



Effect of noise on switches. The stochastic and probabilistic 
analyses, both ultimately based on the chemical master equations 
of the systems 30 , can be compared directly in the third row of 
Figure 1, showing that the stochastic traces at a certain system size 
are compatible with the probability distributions at a smaller system 
size (the relative time scaling is based on the computational 
complexity of the algorithms, and is described in Supplementary 
Materials). This comparison shows close similarity between AM 
and SC, and good similarity between those two and CC apart for 
the fact that CC does not switch fully to x due to the bias provided by 
s. Noticeably, even though a system starts in neutral initial conditions 
(x = b = y), it randomly but rapidly converges to x or y dominance in 
AM, SC, and CC. DC is however quite different from the others, with 
two steep hills representing the stable states, and a wide flat region in 
between where stochastic trajectories may wander. 

The deterministic analysis in the second row or Figure 1, based on 
the mass action ODEs, is not directly comparable with the other two. 
For example, in the extreme case of initial conditions x = y the 
stochastic system still converges with high probability according to 
the 0(N log N) theoretical bound, but the ODEs describe the system 
as remaining in the unstable steady state x = y forever. Therefore, in 
the deterministic case we analyze various non-zero initial discrep- 
ancies between x and y, which still show similarities between AM and 
CC: a logarithmic decrease in disparity leads to linear increase in 
timing. (The deterministic version of DC, omitted, is degenerate: its 
stochastic version is a pure bounded random walk.) 

Both the deterministic and stochastic/probabilistic analyses indi- 
cate that CC is only about twice as slow as the 'optimal' AM algo- 
rithm, and has the same speed as SC. Hence there is not a great 
penalty for the more complex CC network. Further speed optimiza- 
tions of CC are examined in Discussion. 

External control of the AM and CC systems. We now compare the 
steady state response of the four switches to external stimuli. The 
switches of figure 1 are augmented in Figure 2 with external inputs sx 
(switch-to -x) and sy (switch-to-y) acting on both internal steps of the 
pinched arrows in panels B-D. The plots in the second row of figure 2 
integrate three types of steady-state analysis: deterministic, stoch- 
astic, and probabilistic, for different system sizes. In each case we 
vary the input value sx and we observe the steady state output value of 
x, with sy remaining fixed; all rates are 1 .0 and all species start at equal 
level, except for the varying sx and the fixed sy. 

DC shows a hyperbolic response, with no ultrasensitivity, bistabil- 
ity, or hysteresis (Fig. 2A.2). AM and SC show very similar hysteresis 
responses (Fig. 2B,C.2): starting at the origin with increasing sx the 
system stays on the lower trajectory until the jump point indicated by 
the stochastic trace, and continues on the upper trajectory. With 
decreasing sx, the system jumps from high to low at a different jump 
point. CC shows a similar hysteresis response, again somewhat 
depressed in its maxima as in figure 1. The presence of a clear hys- 
teresis cycle in AM, SC, and CC confirms the ability of these systems 
to act as good switches, unlike DC. Although the hysteresis para- 
meters are tunable in the different systems, and AM has fewer species 
and parameters, there is good agreement between their dynamical 
behavior. There is also uniform agreement about the response of each 
network at different systems sizes, with the deterministic, stochastic, 
and probabilistic analysis largely overlapping. Still, the stochastic and 
probabilistic analysis show that noise can induce the transitions 
slightly before the deterministic bifurcation point is reached, thus 
noise can advance the transitions between the two states 31 . 

All the models of figure 1 contain two positive feedback loops. In 
general, it is known that to establish bistability and hysteresis one 
positive feedback (and the presence of nonlinearity) is enough 14 ' 32 . 
Models of the cell-cycle switch have been investigated in face of the 
removal of one of the two positive feedback loops, and it has been 
shown that the removal of one of the two loops decreases the 
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efficiency and the robustness of the CC switch 6,7,12 . In the 
Supplementary Materials we repeat this analysis for the AM and 
CC models and show that the removal of either loops leads to a 
reduced parameter regime for bistability in both models. Thus, also 
in this respect the AM and CC models behave qualitatively similarly, 
but it is worth noting that the CC module is more sensitive (works in 
a smaller parameter range and does not provide full conversion to 
both states), suggesting that AM is indeed the minimal model that 
can provide robust and complete switching dynamics. 

AM and CC switches in the context of oscillators. Having 
established that AM and CC have similar characteristics as stand- 
alone computational systems and dynamical systems, we investigate 
whether AM switches can drive oscillations with characteristics 
similar to those of the cell cycle. There is a large literature on 
various models of cell cycle oscillators 9,10-33 . In many of these a CC 
switch is coupled with another switch in a negative feedback loop 34 . 
Negative feedback loops can induce oscillations even without the 
presence of positive feedbacks, but the positive loops make the 
oscillations more robust as the system periodically switches be- 
tween two states (relaxation oscillations) 35 37 . Connectivity between 
these switches may vary, but all such switches show nonlinearity 38 . 

We consider first an oscillator network structure inspired by a 
mechanical oscillator (the Trammel of Archimedes 39 ), where each 
stable state of a switch produces a transition in the other switch, in a 



repeating pattern (Fig 3A.1 and Supplementary Materials). Two AM 
switches coupled in such way produce a robust limit-cycle oscillator 
(2AM Full) with a self- stabilizing amplitude and frequency 
(Fig. 3A.2). Oscillations are found when ratio between the internal 
AM rates (r ; ) and the external rates connecting the switches (r e ) is 
approximately between 0.2 and 1.0. The detailed response to these 
parameters is investigated in a bifurcation diagram (Fig. 3A.3) and 
on a heat map of the full probability distribution (Fig. 3A.4). We 
observe that even at the trivial parameter combination of all internal 
AM rates equal to 1, the system shows oscillations that are robust in a 
wide regime of external coupling rates that control the strength of the 
negative feedback loops between the two switches. The oscillations 
are born with high amplitude at a SNIC/SNIPER bifurcation 
(Fig. 3A.3), as it has been shown for detailed cell cycle models 11 . 

The model of Figure 3 A contains several feedback loops connect- 
ing the switches, while the majority of cell cycle models contain a 
single negative feedback loop. To better reflect those models, two 
connections have been cut in Figure 3B (2AM Inner) and replaced 
by fixed inputs resembling the effects of cyclin production (sx) and a 
threshold on Cdc20/APC activation (sy) 40 . This network is still able 
to produce robust oscillations, although in a slightly reduced range of 
the parameters of the negative feedback loop (Fig. 3B.3). Some of the 
parameter regimes produce the profiles of relaxation oscillators 
(Fig. 3B.2, top and Fig. 3B.4), typically shown by complex cell-cycle 
models 15,40 . 




Figure 3 | Switches in the context of oscillators. We study the three oscillators depicted in column 1: (A) one consisting of two fully connected AM 
switches; (B) one where two of the connections are replaced by fixed biases; (C) one where we further replace an AM switch with a CC switch. Column 2 
(stochastic high-molecular-count simulations) shows oscillations over time for chosen values of rjr„ where r c is the common rate for the connections 
between switches (gray lines), while r, is the common rate for all the reactions within a switch (black lines). Column 3 explores deterministic parameter 
variation in phase space (xj vs. x 2 ) and bifurcations (r e /r,- or sx vs. Xj). Column 4 combines a deterministic plot (black line), and a probabilistic low- 
molecular-count heat map in phase space (xj vs. x 2 ) for a single value of r e /r;, from column 2. The values of sx and syare fixed to 1/3 of the max value of the 
switching species, and their reactions have rate r e . 
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Finally, in Figure 3C we replace one of the AM switches with a CC 
switch, so that the resulting oscillator (CC-AM Inner) as a whole 
more clearly resembles a cell-cycle oscillator. This also shows that 
even in the context of a wider network the AM and CC switches are 
comparable and interchangeable. The choice of the xJx Y ratio now 
becomes more critical, and large oscillation amplitude is more dif- 
ficult to obtain, but the system still show oscillations similar to those 
of the similarly wired 2AM Inner system and to those of earlier cell 
cycle models 15-40 . The bifurcation diagrams of these two systems also 
look similar to each other and to classical published analysis 11,15 : in 
figure 3B.3 and 3C.3 we plot the steady states and amplitudes of 
oscillations depending on the input parameter sx. 

The external inputs to the oscillators are analogous to cyclin 
production for sx and to cyclin degradation for sy. In natural net- 
works these inputs feed into the oscillator in more complex config- 
urations than in our models 5,10 ' 33 (Supplementary Materials). Still, 
we have demonstrated that switches with the general characteristics 
of AM easily lead to robust oscillators if wired in an appropriate 
configuration. 

Discussion 

The boundary between dynamical systems and computational 
systems is not clear-cut: when does an 'analog' continuous sys- 
tem give rise to 'digital' information processing, and when does a 
collection of ODEs become a digital circuit? Still, there are pre- 
cise characterizations of computational power, both qualitative 
(e.g., whether a computational system is Turing-powerful or not) 
and quantitative (e.g., whether an algorithm runs in 0(log N) or 
not). One can discuss whether a dynamical system falls in one of 
those classes, and what kind of algorithm it implements when 
seen as a computational system. There is increasing interest in 
understanding what kind of physical and biological systems can 



compute in that sense 41 , and it has been proposed that a compu- 
tational outlook can help in our understanding of biological 

systems 42-47 . 

Here we have investigated a classical cell cycle network in terms of 
the algorithm it implements, and we claim we have gained new 
insights in its structure, function, and performance. We have illu- 
strated how the computationally effective AM network and the bio- 
logically relevant CC network are two ways of achieving 'majority 
switching', yielding two distinct algorithmic implementations of 
essentially the same behavior with the same asymptotic properties. 
Moreover, we have shown that their behavior is sufficiently similar 
that they can be readily exchanged for one another in an oscillator 
context. We now discuss some general characteristics of switching 
networks in the context of AM. 

The effectiveness of both AM and CC comes from the fact that 
they both employ two non-linear positive feedback loops. The non- 
linearity in our systems comes from (a purely catalytic version of) 
multi-step modifications 23 . In the context of cell cycle regulation, 
both transcriptional and post-translational modifications have been 
associated with nonlinearity 38 . In relation to the mitotic-switch regu- 
lating post-translational feedback the role of multisite phosphoryla- 
tion has been carefully investigated 48,49 . It has been established that 
for nonlinear response the phosphate groups have to be added and 
removed independently (distributive multisite phosphorylation), 
while if enzymes can add or remove all groups sequentially during 
a single binding event (processive multisite phosphorylation) then 
the nonlinearity is lost 23,50 . AM works like a distributive system 
(Fig. 4A): the common intermediate can be modified independently 
in both directions. (Adding more than one intermediate step for the 
conversion between x and y does not change the fundamental dyna- 
mics of AM.) If instead we use separate intermediates (Fig. 4B), then 
the system works processively: the separate intermediates can be 
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Figure 4 | Nonlinearity in switching systems. We compare the nonlinear dynamics of the AM switch with some alternatives, via sample individual traces 
obtained form stochastic simulations. We have uniformly chosen rate l.Oforunimolecular reactions (in(C) and (D)) and 0.001 for bimolecular reactions, 
with x = y initial conditions (all other species set at zero), and with 15000 max molecule counts as in figure 1. Panel (A) shows the same AM circuit as in 
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(the x-dimer), hence we plot x+2b, the total amount of x. Panel (D) shows direct mutual competition between x and y regarded as two enzymes regulated 
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for y+ b, and hence the system does not operate as a pair of normal enzymatic reactions where those sums are preserved. 
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further modified in only one direction during the transitions between 
x and y, and the system no longer exhibits an irreversible conver- 
gence to stable states. 

Another way to build nonlinearity into biological systems is to let 
the regulators oligomerize, with only the complex having a catalytic 
activity 38 . Similar changes can be made in AM: we can allow dimer- 
ization of x and_y and have these dimers drive transitions between the 
two stable states (Fig. 4C). Models incorporating these changes allow 
AM-like switching between all-x and all-y states, but in the stable 
states we find a mixture of monomers and dimers of the winning 
species. Thus, such protocol would not work efficiently from a com- 
putational point of view. 

Another option is to consider a direct competition network like 
DC but where the catalytic reactions are Michaelis-Menten enzym- 
atic reactions: these are intrinsically nonlinear reactions and could 
potentially replace the nonlinear multisite modification mechanism 
of AM. Note however that again here we have two separate (sub- 
strate-enzyme complex) intermediaries (Fig. 4D). When enzymatic 
reactions work in saturation we have zero-order ultrasensitivity or 
Goldbeter-Koshland switching 51 : these mechanisms are often used in 
biological network models to capture the dynamics of fast switches 52 . 
The AM-like positive feedback network incorporating enzymatic 
reactions, however, does not switch properly as some x and y always 
accumulate in one of the two intermediates 53,54 . On long time scales, 
the primary species and to a lesser degree the complexes oscillate 
randomly, revealing again that the common 'undecided' intermedi- 
ate of the transition between x and y is an important feature of AM. 

Therefore, other plausible mechanisms for achieving nonlinearity 
in chemical systems do not appear to lead to improvements or alter- 
natives to the AM switching algorithm. This suggests that distribu- 
tive multisite modification is the simplest, most effective way of 
introducing nonlinearity in the system. 

Previous literature has shown AM to be an efficient and robust 
decision algorithm 20 and CC to be a robust biological switch 7,12 . We 
used deterministic, stochastic and probabilistic methods to char- 
acterize the kinetic and dynamic behaviors of the two systems and 
those that stand in the way of transforming one into the other. Cells 
do not use AM in its original form, since multifunctional enzymes 
that change their activity through modification are rare 55 . AM is 
robust to parameter changes, but due to its small size it is fragile to 
removal of its components, and thus evolutionary it cannot be 
stable 56,57 . Probably CC is the closest similar system that evolution 
has found to deal with important cellular decisions. 

Similar network structures may have evolved to deal with other 
biological decisions. The cell cycle switch is probably the most well 
characterized biological switch and it is highly conserved among 
eukaryotes. There are many other biological switches that are less 
well characterized 13,58 , but the existence of two positive feedbacks 



in their structure has been proposed. Such systems include the 
membrane transport regulatory switch, were Rab5 helps its activa- 
tion by activating Rabex5 and by inhibiting Rab7 59,60 . Another 
example is the apoptotic switch, where the caspase CASP3 cleaves 
and by this, activates its activator caspases CASP8 and CASP9 while 
also inhibiting its own inhibitor, XIAP 61 . Various models of the cell 
polarity establishment regulatory Cdc42 activation also incorporate 
double positive feedbacks in the system 62 64 . There could be many 
other biological switches that work with similar, Approximate 
Majority-type switching dynamics, and our other findings on the 
AM behavior may apply to these systems as well. 

The main difference between the behavior of AM and CC is that 
CC cannot fully turn on because external biases (s and t on Fig. ID) 
keep inactivating the system. This leads to a reduced maximum in the 
steady state of the CDK-analog species x (Fig. 2D); thus CC is not as 
effective at switching as AM. But nature has solved this problem: 
CDK, through the Greatwall kinase 65 , also inhibits the phosphatase 
PP2 A, which has the catalytic role of both s and t in our CC model. 
With this extra feedback loop, once CDK is active it suppresses the 
activity of the external inhibitors and allows its own full activation. 

The incorporation of this additional feedback into the CC model 
does in fact enable the full activation of x. The extended circuit (GW) 
is shown in Figure 5A, where 5 and t have been merged into s (ana- 
logous to PP2A), with x now inhibiting s. When x is not present, s can 
self-activate (PP2A inactivates its inhibitor, the Greatwall kinase 29 ) 
from an inactive form m, with the help of a fixed bias v. Even though v 
is still continuously counteracting x, its effect is much weaker than 
the direct effect of s in CC, and as a result x can achieve full activation. 
The resulting kinetics is shown in figure 5B for the deterministic 
system (cf. Fig. 2D. 2), and in figure 5C for the stochastic and prob- 
abilistic systems (cf. Fig. 2D. 3). 

Full activation is now achieved, but even more remarkably, the 
switching speed improves in the GW circuit, making it almost as fast 
as AM. This is shown in figure 5D via sample stochastic traces, with 
GW traces coming very close in performance to AM traces as 
opposed to CC traces. The additional circuitry has therefore two 
beneficial effects: improving the activation level of x, and also 
improving the activation speed of x by about a factor of 2, making 
it about as fast as the best known circuit. One could suspect that the 
original decrease in switching speed in SC and CC with respect to 
AM (Fig. 1) is due to the introduction of the intermediate catalysts, 
and therefore that the increased speed in GW is somehow due to the 
direct self-activation of s. However, introducing an intermediate 
species in the s loop still leads to the same speedup (Supplemen- 
tary Materials). 

In conclusion, we have shown that the AM algorithm is the sim- 
plest expression of the task that is performed by the CC network. 
Although the AM algorithm has features that may be biologically 




Figure 5 | Cell cycle switch with Greatwall loop. We study the GW network: a version of CC from Figure 1D.1 enriched with a feedback loop where x 
modulates the sand t biases (now unified into the s species). Panel (B) is analogous to Fig 1D.2, and panel (C) is analogous to Fig 1D.4, both showing an 
improved activation of x. Panel (D) is a comparison between the switching speeds of AM, GW, and CC, showing GW performing better than CC and 
about as well as AM. 
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unrealistic, it can guide the understanding of the properties of more 
complex biological networks. One may also speculate about whether 
related networks, such as SC might exist in nature, and whether 
other complex biological networks may have simple algorithmic 
explanations. 

Methods 

Asymptotic notation. "f{x) is 0(g{x))", or "f is bounded by g asymptotically" is 
defined as: A function over the reals /(x) belongs to the set 0(g{x)) iff there exist real 
numbers k > 0 and x 0 such that/(x) < k ' g(x) for all x > x 0 . Further, "f(x) is co(g(x))", 
or "/dominates g asymptotically" is defined as: A function over the reals f{x) belongs 
to the set oj(g(x)) iff for all real numbers k > 0 there exists a real number x 0 such that 
k • g{x) ^f(x) for all x > x 0 . Note that in, e.g., 0(log x) the basis of the logarithm is not 
important because of the multiplicative constant k and the fact that logarithms in 
different bases differ only by a constant factor. 

Ordinary differential equation models. Implementation of deterministic models 
were done as .ODE models, executable by the XPPAUT (http://www.math.pitt.edu/ 
~- bard/xpp/xpp.html) 26 or Oscill8 (http://sourceforge.net/projects/oscill8) software 
packages. The model files and used parameter sets are provided in Supplementary 
Materials. 

Stochastic models. The various chemical networks were modeled in SPiM Player 
vl.13 http://research.microsoft.com/en-us/projects/spim27. SPiM can represent a 
variety of reactive systems in its modeling language, including chemical reaction 
networks, and can then provide numerical simulations with a Gillespie -style solver: 
the models are included in Supplementary Materials. SPiM's modeling language is a 
general computational language, and it allows us for example to program the slow 
sweeping of sx values in figure 1 within the model and in a single run. In SPiM we can 
study stochastic systems of the order of hundreds of thousands of molecules, and 
obtain an indication of the noise in the systems. For figure 1 row 3 we produced a 
sample of individual Gillespie trajectories: since we used a high molecular count, the 
stochastic noise is quite limited on the individual trajectories (for AM, SC, CC), but 
still the trajectories differ considerably. For figure 2 row 2 we run a Gillespie 
simulation at steady-state for each value of sx, except that all these simulations were 
done in a single run, slowly incrementing sx from 0 to max value, and then slowly 
decrementing it back to 0. The decreasing trajectory was then reflected and 
superimposed on the increasing trajectory, obtaining an image of the hysteresis cycle. 
We used low molecular counts to emphasize the noise. For figure 3 column 2 we run 
stochastic simulations at high molecular count; these are therefore comparable to 
deterministic simulations, but show that even stochastically we can obtain regular 
oscillations. 

Probabilistic models. The various chemical reaction networks were modeled in 
PRISM v4.0.3 http://www.prismmodelchecker.org 28 . PRISM can encode a variety of 
stochastic and probabilistic systems in its modeling language, including stochastic 
chemical networks: the models are included in Supplementary Materials. When 
PRISM is given such a model with initial conditions (initial number of molecules and 
stochastic rates), it generates a continuous -time Markov chain for the system as a 
sparse matrix, where a state of the Markov chain is the vector of number of molecules 
of each species, and a transition in the Markov chain is the propensity of moving from 
the source state to the target state. We can then use the PRISM modelchecking 
language to ask questions about the (possibly very large) Markov chain, such as 
computing the probability of reaching a state from another state. We ask three basic 
questions, and display the answers as heat maps. For figure 1 row 3, what is the 
probability that at time T there will be i molecule of species x ("P = ? [F[T,T]x=i]"); 
for figure 2 row 2, what is the probability that at steady state there will be i molecules of 
species x ("S = ? [x=i]") for each value of sx; and for figure 3 column 4, what is the 
probability that at steady state there will be i molecules of species x 1 and; molecules of 
species x 2 ("S=? [ xl=i &x2=j ]"). The results we can obtain are limited by the size of 
the generated Markov chain; for example, the analysis of the CC circuit in figure 2 
(with 15 molecules of x+b+y, etc.) generates —2.5 million states and —25 million 
transitions. That Markov chain is built in a few seconds, but uses nearly all the 
memory we have available; the computation of the probabilities then takes about 
48 hours on a standard laptop. In contrast, the AM circuit takes only a few minutes 
to analyze completely for the same size. The heat maps were generated in SigmaPlot 
www.sigmaplot.com from the PRISM data. 
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